PPI Example

Label Volume by Trial Type

df_psy <- read.csv("~/Box/my_mri/behavioral/nback_3tb4188_2017_Sep_13_1341_soa.csv") %>%
  rename(onset = start, trial_type = block) %>%
  mutate(duration = ifelse(trial_type == "0-back", 10 * 2.5, 20 * 2.5)) %>%
  group_by(run) %>%
  nest() %>%
  mutate(data = future_map(data, function(x) x %>% select(onset, duration, trial_type)),
         data_trial_type = future_pmap(list(data, tr, n_volumes),  create_trial_type_by_volume_list)) %>%
  unnest(data_trial_type) %>%
  select(-data) %>%
  ungroup()
head(df_psy)

## # A tibble: 6 x 3
##     run volume trial_type
##   <int>  <int> <fct>     
## 1     1      1 fixation  
## 2     1      2 fixation  
## 3     1      3 fixation  
## 4     1      4 3-back    
## 5     1      5 3-back    
## 6     1      6 3-back

Contrast Code

func_fig_contrast <- function(data) {
  data %>%
    as.data.frame() %>%
    mutate(volume = row_number()) %>%
    gather(., psy_contrast, value, -volume) %>%
    ggplot(., aes(volume, value)) +
    geom_line() +
    facet_wrap(~ psy_contrast) +
    theme_minimal()
}

func_add_fig_title <- function(fig, text) {
  fig + 
    labs(title = text)
}

df_psy_contrast <- contrast_code_categorical_variable(df_psy, cbind(nback_vs_baseline = c(1, 1, 1, 1, -4)/5,
                                              task_vs_control = c(-3, 1, 1, 1, 0)/4,
                                              linear_task = c(0, -1, 0, 1, 0),
                                              quadratic_task = c(0, -1, 2, -1, 0)/3)) %>%
  group_by(run) %>%
  nest() %>%
  mutate(data = future_map(data, function(x) x %>% select(contains("psy"))),
         fig = future_map(data, func_fig_contrast),
         fig = future_map2(fig, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_psy_contrast$fig, ncol = 1)

Upsample

df_psy_upsample <- df_psy_contrast %>%
  mutate(data_upsample = future_map(data, function(x) apply(x, 2, function(x) upsample(x, upsample_factor))),
         fig_upsample = future_map(data_upsample, func_fig_contrast),
         fig_upsample = future_map2(fig_upsample, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_psy_upsample$fig_upsample, ncol = 1)

Convolve

df_psy_convolve <- df_psy_upsample %>%
  mutate(data_convolve = future_map(data_upsample, function(x) apply(x, 2, function(x) convolve_afni(x, hrf, tr, n_volumes, upsample_factor))),
         fig_convolve = future_map(data_convolve, func_fig_contrast),
         fig_convolve = future_map2(fig_convolve, paste0("Run:", run), func_add_fig_title))

## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160

ggarrange(plotlist = df_psy_convolve$fig_convolve, ncol = 1)

Downsample

df_psy_downsample <- df_psy_convolve %>%
  mutate(data_downsample = future_map(data_convolve, function(x) apply(x, 2, function(x) downsample(x, upsample_factor))),
         fig_downsample = future_map(data_downsample, func_fig_contrast),
         fig_downsample = future_map2(fig_downsample, paste0("Run:", run), func_add_fig_title))


ggarrange(plotlist = df_psy_downsample$fig_downsample, ncol = 1)

Create Sphere(s)

This step can be skipped if have you have your own mask.

Extract ROI

func_fig_phys <- function(data) {
  colnames(data) <- "data"
  
  data <- as.data.frame(data)
  
  fig <- data %>%
    mutate(volume = row_number()) %>%
    ggplot(., aes(volume, data)) +
    geom_line() +
    theme_minimal()
  return(fig)
}

dir_phys <- "/Volumes/shared/KK_KR_JLBS/Wave1/MRI/FMRI/PPI/Nback_individual_covariates/PHYS_MNI_42_-42_42/3tb1780/"
files_phys <- list.files(dir_phys, "_mean.1D")
path_phys <- paste0(dir_phys, files_phys)
df_phys <- tibble(file = files_phys,
                  path = path_phys) %>%
  mutate(run = str_remove(file, "_mean.1D"),
         run = str_remove(run, "r"),
         data = future_map(path, function(x) read.csv(x, header = F, col.names = "data")),
         fig = future_map(data, func_fig_phys),
         fig = future_map2(fig, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_phys$fig, ncol = 1)

Detrend

df_phys_detrend <- df_phys %>%
  mutate(data_detrend = future_map(data, function(x) detrend(x, 2)),
         fig_detrend = future_map(data_detrend, func_fig_phys),
         fig_detrend = future_map2(fig_detrend, paste0("Run: ", run), func_add_fig_title))


ggarrange(plotlist = df_phys_detrend$fig_detrend, ncol = 1)

Upsample

df_phys_upsample <- df_phys_detrend %>%
  mutate(data_upsample = future_map(data_detrend, function(x) apply(x, 2, function(x) upsample(x, upsample_factor))),
         fig_upsample = future_map(data_upsample, func_fig_phys),
         fig_upsample = future_map2(fig_upsample, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_phys_upsample$fig_upsample, ncol = 1)

Deconvolve

df_phys_deconvolve <- df_phys_upsample %>%
  mutate(data_deconvolve = future_map(data_upsample, function(x) apply(x, 2, function(x) deconvolve_afni(x, hrf))),
         fig_deconvolve = future_map(data_deconvolve, func_fig_phys),
         fig_deconvolve = future_map2(fig_deconvolve, paste0("Run:", run), func_add_fig_title))

## 3dTfitter -RHS /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//data.1D -FALTUNG /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//hrf.1D /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//data_deconvolved.1D 012 -2 -l2lasso -6 
## 3dTfitter -RHS /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//data.1D -FALTUNG /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//hrf.1D /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//data_deconvolved.1D 012 -2 -l2lasso -6 
## 3dTfitter -RHS /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//data.1D -FALTUNG /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//hrf.1D /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_deconvolve_afni//data_deconvolved.1D 012 -2 -l2lasso -6

ggarrange(plotlist = df_phys_deconvolve$fig_deconvolve, ncol = 1)

Multiply/Create Interaction

df_ppi <- tibble(run = c(1:3),
                 data = NA)

for (i in 1:nrow(df_ppi)) {
  df_ppi$data[i] <- list(apply(df_psy_upsample$data_upsample[[i]], 2, function(x) x * df_phys_deconvolve$data_deconvolve[[i]]))
}

df_ppi <- df_ppi %>%
  mutate(fig = future_map(data, func_fig_contrast),
         fig = future_map2(fig, paste0("Run:", run), func_add_fig_title))

ggarrange(plotlist = df_ppi$fig, ncol = 1)

Convolve

df_ppi_convolve <- df_ppi %>%
  mutate(data_convolve = future_map(data, function(x) apply(x, 2, function(x) convolve_afni(x, hrf, tr, n_volumes, upsample_factor))),
         fig_convolve = future_map(data_convolve, func_fig_contrast),
         fig_convolve = future_map2(fig_convolve, paste0("Run:", run), func_add_fig_title))

## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160 
## waver -FILE 0.09375 /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//hrf.1D -input /var/folders/_2/600bt_js7d5dpsn7nbfcpy4x88r01g/T//RtmpaLImsB/temp_convolution_afni//data.1D -numout 4160

ggarrange(plotlist = df_ppi_convolve$fig_convolve, ncol = 1)

Downsample

df_ppi_downsample <- df_ppi_convolve %>%
  mutate(data_downsample = future_map(data_convolve, function(x) apply(x, 2, function(x) downsample(x, upsample_factor))),
         fig_downsample = future_map(data_downsample, func_fig_contrast),
         fig_downsample = future_map2(fig_downsample, paste0("Run:", run), func_add_fig_title))

ggarrange(plotlist = df_ppi_downsample$fig_downsample, ncol = 1)

Initialize

df_design_mat <- tibble(run = c(1:3),
                        data = NA)

for (i in 1:nrow(df_design_mat)) {
  
  temp_psy <- df_psy_downsample$data_downsample[[i]]
  temp_phys <- df_phys_detrend$data_detrend[[i]] %>% rename(phys = residual)
  temp_ppi <- df_ppi_downsample$data_downsample[[i]]
  colnames(temp_ppi) <- str_replace(colnames(temp_ppi), "psy_", "ppi_")

  df_design_mat$data[i] <- list(cbind(temp_psy, temp_phys, temp_ppi))
}

df_design_mat <- df_design_mat %>%
  unnest()

## Warning: `cols` is now required.
## Please use `cols = c(data)`

func_fig_data_heat_map <- function(data) {
  colnames(data) <- paste0(str_pad(1:ncol(data), 2, "left", "0"), "_", colnames(data))
  apply(data, 2, scale_min_max) %>%
    as.data.frame() %>%
    mutate(volume = row_number()) %>%
    gather(., variable, value, -volume) %>%
    ggplot(., aes(variable, volume, fill = value)) +
    geom_raster() + 
    scale_fill_distiller(palette = "Greys", direction = 1) + 
    theme_minimal() +
    theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
    labs(x = NULL)
}

func_fig_data_ts_long <- function(data) {
  colnames(data) <- paste0(str_pad(1:ncol(data), 2, "left", "0"), "_", colnames(data))
  data %>%
  mutate(volume = row_number()) %>%
  gather(., variable, value, -volume) %>%
  ggplot(., aes(volume, value)) +
  geom_line() +
  facet_wrap(~ variable, scales = "free_y", ncol = 1) +
  theme_minimal()
}

func_fig_data_heat_map(df_design_mat)

func_fig_data_ts_long(df_design_mat)

Add nuisance variables

Final

df_design_mat_final <- df_design_mat %>%
  group_by(run) %>%
  nest() %>%
  mutate(data = future_map(data, function(x) x %>% mutate(time_linear = row_number(), 
                                                          time_linear = scale_min_max(time_linear),
                                                          time_quadratic = time_linear^2))) %>%
  unnest() %>%
  ungroup()

## Warning: `cols` is now required.
## Please use `cols = c(data)`

func_fig_data_heat_map(df_design_mat_final)

func_fig_data_ts_long(df_design_mat_final)